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Abstract. Semi-discrete Galerkin formulations of transient wave equations, either with con- 
forming or discontinuous Galerkin finite element discretizations, typically lead to large sys- 
tems of ordinary differential equations. When explicit time integration is used, the time-step is 
constrained by the smallest elements in the mesh for numerical stability, possibly a high price 
to pay. To overcome that overly restrictive stability constraint on the time-step, yet without 
resorting to implicit methods, explicit local time-stepping schemes (LTS) are presented here 
for transient wave equations either with or without damping. In the undamped case, leap-frog 
based LTS methods lead to high-order explicit LTS schemes, which conserve the energy. In 
the damped case, when energy is no longer conserved, Adams-Bashforth based LTS methods 
also lead to explicit LTS schemes of arbitrarily high accuracy. When combined with a finite 
element discretization in space with an essentially diagonal mass matrix, the resulting time- 
^-j- marching schemes are fully explicit and thus inherently parallel. Numerical experiments with 
continuous and discontinuous Galerkin finite element discretizations validate the theory and 
illustrate the usefulness of these local time-stepping methods. 
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1 Introduction 



The efficient numerical simulation of transient wave phenomena is of fundamental im- 
portance in a wide range of applications from acoustics, electromagnetics or elasticity. 
Although classical finite difference methods remain a viable approach in rectangu- 
lar geometry on Cartesian meshes, their usefulness is quite limited in the presence of 
complex geometry, such as cracks, sharp corners or irregular material interfaces. In 
contrast, finite element methods (FEMs) easily handle unstructured meshes and local 
refinement. Moreover, their extension to high order is straightforward, a key feature to 
keep numerical dispersion minimal. 

Semi-discrete finite element Galerkin approximations typically lead to a system of 
ordinary differential equations. However, if explicit time-stepping is subsequently em- 
ployed, the mass matrix arising from the spatial discretization by standard conforming 
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finite elements must be inverted at each time-step: a major drawback in terms of effi- 
ciency. To overcome that difficulty, various "mass lumping" techniques have been pro- 
posed, which effectively replace the mass matrix by a diagonal approximation. While 
straightforward for piecewise linear elements @|32), mass lumping techniques require 
special quadrature rules and additional degrees of freedom at higher order to preserve 
the accuracy and guarantee numerical stability lfT0ll20l . 

Discontinuous Galerkin (DG) methods offer an attractive and increasingly popular 
alternative for the spatial discretization of time-dependent hyperbolic problems 0] |7J 
ET]|22l|30l[33. Not only do they accommodate elements of various types and shapes, 
irregular non-matching grids, and even locally varying polynomial order, and hence 
offer greater flexibility in the mesh design. They also lead to a block-diagonal mass 
matrix, with block size equal to the number of degrees of freedom per element; in 
fact, for a judicious choice of (locally orthogonal) shape functions, the mass matrix is 
truly diagonal. Thus, when a spatial DG discretization is combined with explicit time 
integration, the resulting time-marching scheme will be truly explicit and inherently 
parallel. 

In the presence of complex geometry, adaptivity and mesh refinement are certainly 
key for the efficient numerical simulation of wave phenomena. However, locally re- 
fined meshes impose severe stability constraints on explicit time-marching schemes, 
where the maximal time- step allowed by the CFL condition is dictated by the small- 
est elements in the mesh. When mesh refinement is restricted to a small region, the 
use of implicit methods, or a very small time-step in the entire computational domain, 
are a very high price to pay. To overcome this overly restrictive stability constraint, 
various local time-stepping (LTS) schemes lfl2l[T3l[T71l were developed, which use ei- 
ther implicit time-stepping or explicit smaller time-steps, but only where the smallest 
elements in the mesh are located. 

Since DG methods are inherently local, they are particularly well-suited for the de- 
velopment of explicit local time-stepping schemes ll30l . By combining the sympletic 
Stormer-Verlet method with a DG discretization, Piperno derived a symplectic LTS 
scheme for Maxwell's equations in a non-conducting medium ll35l . which is explicit 
and second-order accurate. In ll34ll . Montseny et al. combined a similar recursive 
integrator with discontinuous hexahedral elements. Starting from the so-called arbi- 
trary high-order derivatives (ADER) DG approach, alternative explicit LTS methods 
for Maxwell's equations [39 1 and for elastic wave equations [18| were proposed. In 
|[T9"ll , the LTS approach from Collino et al. |[T2l [131 was combined with a DG-FE 
discretization for the numerical solution of symmetric first-order hyperbolic systems. 
Based on energy conservation, that LTS approach is second-order and explicit inside 
the coarse and the fine mesh; at the interface, however, it nonetheless requires at every 
time-step the solution of a linear system. More recently, Constantinescu and Sandu 
devised multirate explicit methods for hyperbolic conservation laws, which are based 
on both Runge-Kutta and Adams-Bashforfh schemes combined with a finite volume 
discretization lfl4l [T31 . Again these multirate schemes are limited to second-order 



Explicit local time-stepping methods 



3 



accuracy. 

Starting from the standard leap-frog method, Diaz and Grote proposed energy con- 
serving fully explicit LTS integrators of arbitrarily high accuracy for the classical 
wave equation [16]; that approach was extended to Maxwell's equations in 11251 for 
non-conductive media. By blending the leap-frog and the Crank-Nicolson methods, a 
second-order LTS scheme was also derived there for (damped) electromagnetic waves 
in conducting media, yet this approach cannot be readily extended beyond order two. 
To achieve arbitrarily high accuracy in the presence of dissipation, while remaining 
fully explicit, explicit LTS methods for damped wave equations based on Adams- 
Bashforth (AB) multi-step schemes were proposed in Il26ll - see also Il27l . They can 
also be interpreted as particular approximations of exponential- Adams multistep meth- 
ods on . 

The rest of the paper is organized as follows. In Section 2, we first recall the standard 
continuous, the symmetric interior penalty (IP) DG and the nodal DG formulations. 
Next in Section 3, we consider leap-frog based LTS methods, both for the undamped 
and the damped wave equation. In the undamped case, we show how to derive explicit 
LTS methods of arbitrarily high order; these methods also conserve a discrete version 
of the energy. In the damped case, we present a second-order LTS method by blending 
the leap-frog and the Crank-Nicolson scheme; however, this approach does not eas- 
ily extend to higher order. To achieve arbitrarily high accuracy even in the presence 
of dissipation, we then consider LTS methods based on Adams-Bashforth multi-step 
schemes in Section 4. Finally in Section 5, we present numerical experiments in one 
and two space dimensions, which validate the theory and underpin both the stability 
properties and the usefulness of these high-order explicit LTS schemes. 

2 Finite element discretizations for the wave equation 

We consider the damped wave equation 

u tt + out - V • (c 2 V«) = / in D. x (0, T) , 

u(-,t) =0 on dQ. x (0,T), (2.1) 
u(-,0) = uq , Ut(-,0) = v inQ., 

where Q. is a bounded Lipschitz domain in R d , d = 1 , 2, 3. Here, / G L 2 (0, T; L 2 (Q.)) 
is a (known) source term, while uq 6 Hq(Q) and vq 6 L 2 (D.) are prescribed initial 
conditions. At the boundary, d£l, we impose a homogeneous Dirichlet boundary con- 
dition, for simplicity. We assume that the damping coefficient, a = a(x) and the speed 
of propagation c = c(x) are piecewise smooth and satisfy the bounds 

< <j(x) < a* < oo , < < c(x) < c* < oo , x 6 £1 ■ 



We shall now discretize ( 2. 1 1 in space by using any one of the following three dis- 
tinct FE discretizations: continuous (H l -conforming) finite elements with mass lump- 
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ing, a symmetric IP-DG discretization, or a nodal DG method. Thus, we consider 
shape -regular meshes Th that partition the domain D, into disjoint elements K, such 
that D, = UxeT h K- The elements are triangles or quadrilaterals in two space dimen- 
sions, and tetrahedra or hexahedra in three dimensions, respectively. The diameter of 
element K is denoted by hx and the mesh size, h, is given by h = max^- e f h hpc- 



2.1 Continuous Galerkin formulation 

The continuous (H l -conforming) Galerkin formulation of (2.1 1 starts from its weak 
formulation: find u e [0, T] -> Hq (Q) such that 

(uu,tp) + ((TUt,<p) + (cVu,cVtp) = (f,cp) V^effoHn), te(0,T), 

(2. 2 j 

u(-,0) =u Q , ut(-,0) =v , 



where (■, ■) denotes the standard L 2 inner product over D,. It is well-known that (2.2 1 
is well-posed and has a unique solution Il33l . 

For a given partition Th of £2, assumed polygonal (2d) or polyhedral (3d) for sim- 
plicity, and an approximation order £ > 1, we shall approximate the solution u(-, t) of 
(2.2 1 in the finite element space 

:= e ffi(Q) : p\ K £ S\K) VKeTh) , 

where <S (-fT) is the space V (K) of polynomials of total degree at most £ on K if if 
is a triangle or a tetrahedra, or the space Q e (K) of polynomials of degree at most £ in 
each variable on if if if is a parallelogram or a par allele piped. Here, we consider the 
following semi-discrete Galerkin approximation of (2.2 1: find u h : [0,T] — > V h such 
that 

(4 t ,v) + (<Ju>?,<p) + (cX7u h ,cX7ip) = (f,<p) Vip£V h , te(0,T), 

^(■,0) = n ftUo , ^(-,o) = n h «o. 

Here, denotes the L 2 -projection onto V h . 

The semi-discrete formulation (2.3 1 is equivalent to the second-order system of or- 
dinary differential equations 

d 2 XJ 

M ¥ (t) + Mff ¥ (t) + KU(t)=F(t) ' te (0 ' T) ' 

MU(0)=«§, M^(0) = « \ 

Here, U denotes the vector whose components are the time-dependent coefficients of 
the representation of u h with respect to the finite element nodal basis of Vh, M the 
mass matrix, K the stiffness matrix, whereas denotes the mass matrix with weight 
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a. The matrix M is sparse, symmetric and positive definite, whereas the matrices K 
and M CT are sparse, symmetric and, in general, only positive semi-definite. In fact, K 



is positive definite, unless Neumann boundary conditions would be imposed in ( |2.1[ > 
instead. Since we shall never need to invert K, our derivation also applies to the semi- 
definite case with purely Neumann boundary conditions. 

Usually, the mass matrix M is not diagonal, yet needs to be inverted at every time- 
step of any explicit time integration scheme. To overcome this diffculty, various mass 
lumping techniques have been developed iT0l[TTl l8ll9l. which essentially replace M 
with a diagonal approximation by computing the required integrals over each element 
K with judicious quadrature rules that do not effect the spatial accuracy (4]|. 

2.2 Interior penalty discontinuous Galerkin formulation 



Following [21] we briefly recall the symmetric interior penalty (IP) DG formulation 
of ( |2.1[ >. For simplicity, we assume in this section that the elements are triangles or 
parallelograms in two space dimensions and tetrahedra or parallelepipeds in three di- 
mensions, respectively. Generally, we allow for irregular (fc-irregular) meshes with 
hanging nodes [5|. We denote by £^ the set of all interior edges of Th, by £j? the set 
of all boundary edges of Th, and set £h = Eh U £f L ' . Here, we genetically refer to any 
element of Eh as an "edge", that is a real edge in 2d and a face in 3d. 

For a piecewise smooth function <p, we introduce the following trace operators. Let 
e 6 £? be an interior edge shared by two elements K + and K~ with unit outward 
normal vectors n ± , respectively. Denoting by v ± the trace of v on dK taken from 
within K , we define the jump and the average on e by 

[<p\ := <p + n + + <p-zT , l<p} := {<p + + <p~)/2. 

On every boundary edge e 6 £&, we set [yj] := i^n and ^ipj := tp. Here, n is the 
outward unit normal to the domain boundary d£l. 

For a piecewise smooth vector- valued function ip, we analogously define the average 
across interior faces by ^ipj ■= + tp~)/2, and on boundary faces we set ■jf?/'}} : = 
ip. The jump of a vector-valued function will not be used. For a vector-valued function 
ip with continuous normal components across a face e G Eh, the trace identity 



(p + (n + ■ + <p (n -if) 



on e . 



immediately follows from the above definitions. 

For a given partition Th of £2 and an approximation order I > 1 , we wish to approx- 



imate the solution u(t, ■) of (2.1 1 in the finite element space 



V h := e L 2 (Q) : V \ K e S e (K) VK G 77,} , 

where S e (K) ist the space V (K ) (for triangles or tetrahedra) or Qr(K ) (for quadrilat- 
erals or hexahedra). Thus, we consider the following (semidiscrete) DG approximation 
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of (2.1 1: find u h : [0, T] -> V h such that 



(4, y> + (<n#, V ) + a ft (ti h , <p) = tf,<p) V<peV h , te (0, T , 

u h (-,o) = n h uo, u?(-,o) = n^ . 

Here, II/j again denotes the L 2 -projection onto I/' 1 whereas the DG bilinear form 
a h (-, ■), defined on V h x V h , is given by 



V?) := 53 / c 2 Vu-Vipdx-J2 I M ■ {{c 2 v ^}} d ^ 
J] /M-^V^^+J] / a[u] • [<p]dA. 



(2.6) 



The last three terms in (2.6 1 correspond to jump and flux terms at element boundaries; 



they vanish when u,<p, 6 Hq(Q.) n H 1+m (D.) for m > i. Hence, the above semi- 



discrete DG formulation (|2.5[) is consistent with the original continuous problem ( 2.2 1 



In (2.6 1 the function a penalizes the jumps of u and v over the faces of 7~h- To define 



it, we first introduce the functions h and c by 

I min{h K i , h K -}, e e , \ max{c\K+(x),c\ K -(x)}, ee^ 



h K , e64 B , c\k{x), eeS h 



Then, on each e 6 Ey t , we set 



acV, (2.7) 



where a is a positive parameter independent of the local mesh sizes and the coefficient 
c. There exists a threshold value a m i n > 0, which depends only on the shape regularity 
of the mesh and the approximation order i such that for a. > a m i n the DG bilinear 
form ah is coercive and, hence, the discretization is stable [121 El . Throughout the rest 
of the paper we shall assume that a > a m i n so that the semi-discrete problem ( |2.5[ > 
has a unique solution which converges with optimal order Il2ni22ll23ll24l . In ll2~Tll2"4ll . 
a detailed convergence analysis and numerical study of the IP-DG method for ( |2.6[ > 
with (7 = was presented. In particular, optimal a-priori estimates in a DG-energy 
norm and the L 2 -norm were derived. This theory immediately generalizes to the case 



a > 0. For sufficiently smooth solutions, the IP-DG method (2.6 1 thus yields the 
optimal L 2 -error estimate of order 0(h l+x ). 

The semi-discrete IP-DG formulation (|2.5| is equivalent to the second-order system 



of ordinary differential equations ( 2.4 1. Again, the mass matrix M is sparse, symmetric 
and positive definite. Yet because individual elements decouple, M (and M<j) is block- 
diagonal with block size equal to the number of degrees of freedom per element. Thus, 
M can be inverted at very low computational cost. In fact, for a judicious choice of 
(locally orthogonal) shape functions, M is truly diagonal. 
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2.3 Nodal discontinuous Galerkin formulation 

Finally, we briefly recall the nodal discontinuous Galerkin formulation from QUI for 
the spatial discretization of (2.1 1 rewritten as a first-order system. To do so, we first let 
v := uu w := -V«, and thus we rewrite (2.1 1 as the first-order hyperbolic system: 

v t + crv + V ■ (c 2 w) = / in Q x (0, T) , 

Wi + X7v = in D. x (0, T) , 

v(-,t) = 0, w(-,t) = on 9(2x(0,T), 

«(-,0) = vo, w(-,0) = -Vtto in ^> 
or in more compact notation as 

qi + Sq + V-7(q) = S, 

with 



(2.8) 



(2.9) 



q 



2 T 
C W 



Vldxd 





Following [ 30 1, we now consider the following nodal DG formulation of (2.9 1: find 
q h : [0, T] -> V fe such that 

(qf^) + (Sq h ,^)+^(q^,V) = (S,V) V^eV ft , te(0,T). (2.10) 
Here denotes the finite element space 

V h := G L 2 (£2) d+1 : G ViT G 

for a given partition 7ft of D, and an approximation order £ > 1 . The nodal-DG bilinear 
form aft(-, •) is defined on V™ x V' 1 as 



aft(q,^) 



^ / (V--F(q))-^- £ [ (n-T(<D-(n-T(q))*)-ipdA, 

Ken Jk eee h Je 



where (n ■ .F(q))* is a suitably chosen numerical flux in the unit normal direction n. 
The semi-discrete problem (2.10l has a unique solution, which converges with optimal 
order in the L 2 -norm ll30l . 

The semi-discrete nodal DG formulation ( 2. 10 1 is equivalent to the first-order sys- 
tem of ordinary differential equations 



U + M a Q(t) + C Q(t) = F(t) 



te(o,r). 



(2.11) 



Here Q denotes the vector whose components are the coefficients of q h with respect to 
the finite element basis of V" and C the DG stiffness matrix. Because the individual 
elements decouple, the mass matrices M and M CT are sparse, symmetric, positive semi- 
definite and block-diagonal. Moreover, M is positive definite and can be inverted at 
very low computational cost. 
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3 Leap-frog based LTS methods 

Starting from the well-known second-order "leap-frog" scheme, we now derive an ex- 
plicit second-order LTS scheme for undamped waves. By using the modified equation 
approach, we then derive an explicit fourth-order LTS method for undamped waves. 
Finally, by blending the leap-frog and the Crank-Nicolson methods, we also present a 
second-order LTS scheme for damped waves. 
We consider the semi-discrete model equation 



M^(t) + KU = F(t), 
at 



(3.1) 



where M and M a are symmetric positive definite matrices and K is a symmetric pos- 
itive semi-definite matrix. Moreover, we assume that the mass matrix M is (block-) 
diagonal, as in ( |2.4) . We remark, however, that the time integration techniques pre- 
sented below are also applicable to other spatial discretizations of the damped wave 
equation that lead to the same semi-discrete form p. I} . 

Because M is assumed essentially d iago nal, can be explicitly computed and 
inverted at low cost. Thus, we multiply (3.1 1 by M~2 to obtain 



d 2 z , s ^ dz , . 



Az(t) = R(i) , 



(3.2) 



with z = MsU, D = M"2M g M"5, A = M"2KM"2 and R = M"5F. Note 
that A is also sparse and symmetric positive semidefinite. For undamped waves, D 
vanishes and hence energy is conserved, whereas for damped waves D is nonzero and 
energy is dissipated. We shall distinguish these two situations in the derivation of local 
time-stepping schemes below. 



3.1 Second-order method for undamped waves 



For undamped waves, ( |3.2| reduces to 

d?z 
dt 2 



Az = R. 



(3.3) 



Since for any / £ C 2 , we have 

f(t + At)-2 f(t) + f(t - At) = At 2 J (1 - \e\)f"(t + 9 At) dd , (3.4) 



the exact solution z(t) of (3.3 1 satisfies 



z{t + At) -2z(t) +z(t- At) = At 2 J (1 - |0|) (R(t + 0At) - Az{t + 6At)) 



d9. 
(3.5) 
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The integral on the right side of (3.5 1 represents a weighted average of R(s) — A z(s) 
over the interval [t — At, t + At], which needs to be approximated in any numerical 
algorithm. If we approximate in (3.5 1 Az(t + 8 At) and R(t + 9 At) by Az(t) and 
R(f), respectively, and evaluate the remaining ^-dependent integral, we obtain the 
well-known second-order leap-frog scheme with time-step At, 

Zn+i ~ 2z„ + z n _! = At 2 (R„ - Az n ) , R„ ~ R(t„), z n ~ z(t„) , (3.6) 

which, however, would require At to be comparable in size to the smallest elements in 
the mesh for numerical stability. 

Following Ifl6ll28l . we instead split the vectors z(t) and R(t) as 

z(t) = (I - P)z(t) + Pz(t) = z^'^Ut) + Z [ fine l(t) , 

(3.7) 

R(t) = (I - P)R(t) + PR(t) = R [coarse] (t) + R[ fine l(t) , 

where the projection matrix P is diagonal. Its diagonal entries, equal to zero or 
one, identify the unknowns associated with the locally refined region, where smaller 
time-steps are needed. To circumvent the severe CFL restriction on At in the leap- 
frog scheme, we need to treat z[ flne l(t) and R[ flne l(t) differently from z[ coarse l(t) and 
R[ coarse l(t) in 

z(f + At) - 2 z(t) + z(t - At) 

= At 2 ^ (1 - |0|) {R[ coarse l(t + 9 At) + R[ fine l(t + 9 At) (3 _ 8) 

-A (> oarse ] (t + 9 At) + z [ flne l (t + 9 At))} d8 . 

Since we wish to use the stand ard leap-frog scheme in the coarse part of the mesh, we 
approximate the terms in (^ that involve z [ coarse l (t + 9 At) and R[ coaise l (t + 9 At) by 
their values at t, which yields 

z(t + At) - 2 z(t) + z(t - At) ~ At 2 {(I - P)R(f) - A(I - P)z(t)} 

r l (3-9) 
+ At 2 J (1 - \9\) {PR(t + 9At) - APz(t + 6»At)} d9 . 

Note that A and P do not commute. 

Next for fixed t, let z(r) solve the differential equation 

0(r) = (I - P)R(t) - A(I - P)z(t) + PR(t + r) - APz(r) , ^ 

z(0)=z(t), z'(0) = v, 

where v will be specified below. Again from ( |3.4) , we deduce that 

z(At) - 2 z(0) + z(-At) = At 2 {(I - P)R(t) - A(I - P)z(t)} 

ri (3.11) 
+ At 2 J (1 - \9\) {PR(t + 0At) - APz(6»At)} dS . 
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From the comparison of (3.9 1 with ( |3.11[ >, we infer that 

z(t + At)-2 z{t) + z(t - At) ~ z(At) - 2 z(0) + z(-At) , 
or equivalently 

z(t + At) + z(t - At) ~ z(At) + z(-At) . (3.12) 
In fact from Taylor expansion and ( |3.3| >, we obtain 

z(At) + z(-At) = 2z(0) + z"(0)At 2 + 0(At 4 ) 

= 2z(t) + (R(t) - Az(t))At 2 + 0(At 4 ) = z{t + At) + z(t - At) + 0{At 4 ) . 

Thus to advance z(t) from t to t + At, we shall evaluate z(At) + z(— Ai) by solving 
( 3.10| numerically. 

To take advantage of the inherent symmetry in time and thereby reduce the compu- 
tational effort even further, we now let 

q(r) = z(r) + z(-t) . 

Then, q(r) solves the differential equation 



dr 2 



(r) = 2 {(I - P)R(i) - A(I - P)z(i)} + P {R(t + r) + K{t - r)} 



- APq(r) , 
q(0) =2z(t), q'(0) =0, 



(3.13) 



with 

z(t + At) + z(t - At) = q(Af) + 0(At 4 ) . (3.14) 

Note that q(At) does not depend on the value of v. Now, we shall approximate the right 
side of ( 3.5 • by solving (3.13 • on [0,At], and then use (3.14l to compute z(t + At). 
Thus, we need the numerical value of q(r) only at At. 

In summary, the second-order LTS algorithm for the solution of (3.3 1 computes 
z n+ i ~ z(t + At), given z„ and z„_i, as follows: 

LTS-LF2(jj) Algorithm 

(i) Set w := (I — P)R n — A(I — P)z„ and qo := 2z„. 



(ii) Compute qj/j, := qo 



1 /AA 



(2w + 2PR„, - APqo) . 

P J 



(iii) For m = 1, . . . ,p — 1, compute 

Cl(rre+l)/p : = 2q m / p — q( m -l)/p + 



(2w + P(Rn,m + R-k,- 



APq ra / p ) 
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(iv) Compute z n+ i := -z„_i + qi. 

Here, we have used the notations R„ jm — R(i n + r m ) and R nj _ m ~ R(t n — r m ), 
where t n = nAt and r m = raAr; note that R n .o — R{t n + To) = R(i n ) — R ?l - Steps 
1-3 correspond to the numerical solution of ( |3.13} until r = At with the leap-frog 
scheme, using the local time-step At = At /p. For P = 0, that is without any local 
time-stepping, we recover the standard leap-frog scheme. If the fraction of nonzero 
entries in P is small, the overall cost is dominated by the computation of w, which 
requires one multiplication by A(I — P) per time-step At. All further matrix-vector 
multiplications by AP only affect those unknowns that lie inside the refined region, 
or immediately next to it. 



Proposition 3.1. For R(t) 6 C 2 ([0, T]), the local time-stepping method LTS-LF2(p) 
is second-order accurate. 



Proof. See El. □ 
To establish the stability of the LTS-LF2(p) scheme we consider the homogeneous 



case, R„ = 0. Then, the standard leap-frog scheme ( 3.6 1 conserves the discrete energy 

1 " 



2 



At \ z„ + i — z„ z n+ i — z n \ / z n+ \ + z n z n+ i + z„ 
X A J At ' At / + \ 2 ' 2~ 



(3.15) 

Here E n+ 2 ~ E(t , i ) and the angular brackets denote the standard Euclidean inner 



product. Since A is symmetric, the quadratic form in (3.15 1 is also symmetric. For 
sufficiently small At it is also positive semidefinite and hence yields a true energy. 

To derive a necessary and sufficient condition for the numerical stability of the LTS- 
LF2(p) scheme, we exhibit a conserved discrete energy for the LTS-LF2(p) algorithm 
with R„ = 0. Following lfl6l . we first rewrite the LTS-LF2(p) scheme in "leap-frog 
manner". 



Proposition 3.2. The local time-stepping scheme LTS-LF2(p) with R njm = is equiv- 
alent to 



where A„ is defined by 



z ra+l — 2z„ Z n _i At A.pZ n , 



2 ^/A^ 2j 



3=1 



A, = A-4]T - a*(AP)'A (3.16) 
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and the constants op- are given by 



a\ = \, ot\ = 6, a\ - 

p+l 2 i o p p—i 

a, = m + 2a\ — ct\ , 



v P+1 

j 

p+l 



2a i ~ a i 1 ~ a j-i' J = 2 ' • ■ ' > p ~ 2 ' 



a;_! =2a 3 ,_ 1 -a p _ 2 , 
u p ~ u p-i ■ 
Furthermore, the matrix A p is symmetric. 

Proof. See lfT6l and ll25l. 



Proposition 3.3. The local time-stepping scheme LTS-LF2(p) with R„ = conserves 
the energy 



E" 



1 



4 



At 



At 



z n+l z n z n+l z n \ _^ / z n+l H~ z n z n+l "I - z r, 



(3.17) 



Proof. By symmetry of A p , this standard argument is similar the proof of (3.15 1; see 
also |[T6l for details. □ 



As a consequence, the LTS-LF2(p) is stable if < (At 2 /4)A max (A J ,) < 1; note that 
the matrix A p itself also depends on At. 



3.2 Fourth-order method for undamped waves 

In the absence of damping, the wave equation corresponds to a separable Hamiltonian 
system. This fact explains the success of symplectic integrators, such as the Stormer- 
Verlet or the leap-frog method, when combined with a symmetric discretization in 
space. Indeed the fully discrete numerical scheme will then conserve (a discrete ver- 
sion of) the energy, too. Clearly, standard symplectic partitioned Runge-Kutta (Lobatto 
IIIA-IIIB pairs) or composition methods 11281 can be used to achieve higher accuracy 
ll36ll . Because the Hamiltonian here is separable, those higher order versions will also 
remain explicit in time, like the Stormer-Verlet method. Since damped wave equations 
are linear, we instead opt for the even more efficient modified equation (ME) approach 
ll38ll in this section, which leads to explicit LTS of arbitrarily high (even) order. 



Following the ME approach, we replace Az(t+8 At) in ( 3.5 1 by its Taylor expansion 



/ (fi At 2 (fi At^ \ 

Az(t + 6 At) = A ( z(t) + 6 At z'(t) + 7." (t) + — z"'(t) J + 0(At 4 
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Then, the integrals involving odd powers of 9 vanish. Next, by using that z"(t) = 
R(i) — Az(t) and the Simpson quadrature rule for the term that involves R(f + 8 At), 
we obtain the fourth-order modified equation scheme. 



5 m+l 



2z m , 
At 2 



z m— 1 



R»? 
l 



Az, 



At 2 



A 2 2 



At 1 
12 



AR„ 



(3.18) 



+ ~ (R-m-l/2 ~ 2R m + Rm+1/2) + C(Ai ) , 



where z m ~ z(t m ), R m ~ R(i m ) and R m ±i/2 — P-C^m^ Ai/2). Clearly, integration 
schemes of arbitrary (even) order can be obtained by using additional terms in the Tay- 
lor expansion. Since the maximal time-step allowed by the fourth-order ME method 
is about 70% times larger than that of the leap-frog scheme ifTOl . the additional work 
needed for the improved accuracy is quite small; hence, the ME method is extremely 
efficient. 



We now derive a fourth-order LTS method for (3.3 1. Similarly to the derivation in 
Section 3.1 we split the vectors z(t) and R(£) in (3.8 1 into a fine and a coarse part, 
and shall treat z^(t) and R[ fine l(i) differently from z^ oar ^(t) and R[ coarse l(i). We 
expand z [ coarse l (t + 6 At) in Taylor series as 



\t + 6 At) = z^ Kt \t) + 6 At 



dz^ 



(*) + 



e 2 At 2 d 2 z [ coarse i 



dt 

e 3 At 3 d 3 z[ coarse 



dt 3 



2 dt 2 

(t) + 0(At 4 ) 



(0 



and insert it into (3.8 1. In (3.8 1, the integrals involving odd powers of 6 vanish. By 
using 



d 2 z 



coarse] 

dW 



dh 



(t) = (I-P)^-(t) = (I-P)R(t) 



-P)Az(i) 

and the Simpson quadrature rule for the term in ( |3.8| that involves R[ coarse l(i + OAt), 
we find that 

z(t + At) - 2 z(t) + z(t - At) 



Af 4 

At 2 {(I - P)R(t) - A(I - P)z(t)} + — A(I - P)Az(t) 



At 4 At 2 
-_A(I-P)R(t) + 1 -(I-P)<jR(/ 



At 

y 



2R(t) + R [t 



(3.19) 



+At 2 J 1 (1 - \6\) {R [fine] (i + OAt) - Az[ flne '(i + 6At)} d6 . 



Hence, if P = we recover the standard ME scheme (|3.18|). 
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Similarly to Section 3. 1 we now approximate the right-hand side of ( 3. 19 1 by solv- 
ing the following differential equation for z(r) 



fz 

dr 2 



r) = (I-P)R(i)-A(I-P)z(t) 



: (I-P) R t 



At 
2~ 



2R(t) + R t 



At 
2~ 



2 2 

+ jA(I - P)Az(t) - yA(I - P)R(t) + PR(t + r) - AP5(r) , 
2(0) = z(t), z'(0) =i/, 
where i/ will be specified below. Again, using Taylor expansions, we infer that 
z(t + At) + z(t - At) = z(At) + z(-At) + 0{At 6 ) . 

Again, the quantity z(At) + z(— At) does not depend on the value of v, which we set 
to zero. As in Section [3~Tj we set q(r) := z(r) + z(— r), which solves the differential 
equation 



d^ 1 



r)=2{(I-P)R(t)-A(I-P)z(t)} 



+ -(I-P)<|R(t 



At 

y 



2R(t) + R t 



At 

T 



+ r 2 A(I - P)Az(t) - r 2 A(I - P)R(t) 
+ P{R(t + r)+R(t-r)} 
- APq(r) , 
q(0) =2z(t), q'(0)=0. 



(3.20) 



Thus, we have 



z(t + At) + z(t - At) = q(At) + C(At 6 



(3.21) 



Now, we approximate the right side of <\3.5\ by solving ( |3.20[ ) with the fourth-order 
ME method on [0, At] with a smaller time step At = At/p, and then use (3.21 1 to 
compute z(t + At). 

In summary, the fourth-order LTS algorithm for the solution of (3.3 1 computes 
z n+ i ~ z(t + At), given z„ and z n _i, as follows: 

LTS-LFME4(p) Algorithm 

(i) Set q := 2z„, Wl := (I - P)R„ - A(I - P)z„, 

w 2 ■= A(I - P)Az n - A(I - P)R„ andri := R„_i/ 2 - 2R n + R„ +1/2 - 
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(ii) Compute 



u := 2wi + - (I - P)n + 2PR Jlj0 - APq 

1 fAt\ 2 1 /AA 4 / _ / 2 x : 



q 1/p :=q + ^ 7 j ^^[j 2W2 + 2 UJ Pri " AP N^ 



(iii) For m = 1, . . . ,p — 1, compute 



fmAt\ 



2 



Wl := 2w! + P)n + ( ) w 2 + P (Rn,m + R-™,-m) - APq m/p , 



r :— R„ jm _i/2 — 2R„ iT7l + R njm+ i/2 + R nj _ m _i/2 
2R n __ m + R n; _ m +i/2 ) 

w 2 := 2w 2 + (-^) Pr-APui, 
\mAt / 

At\ 2 1 /At x4 



q(m+i)/ P : = 2qm/ P - q(m-i)/p + ( — ) «i + ^ ( ^ ) M2 • 



(iv) Compute z n+l := — z n _i + qi. 



Here, Steps 1-3 correspond to the numerical solution of (3.20 1 until r = At with 
the ME approach using the local time-step At = At/p. The LTS-LFME4(p) algo- 
rithm requires three - two, without sources - multiplications by A(I — P) and 2p 
further multiplications by AP. For P = 0, that is without any local time-stepping, the 



algorithm reduces to the modified equation scheme (3.18 1 above. 



3.3 Second-order leap-frog/Crank-Nicolson based method for damped 

waves 

We shall now derive a second-order LTS method for (|3.2| in a general form with 



D 7^ 0. In contrast to the time-stepping scheme presented in Section 3.1 for the case 
D = 0, we are now faced with several difficulties due to the additional Dz'(f) term. 
First, we shall treat that term implicitly to avoid any additional CFL restriction; else, 
the stability condition will be more restrictive than that with the LTS-LF2(p) scheme, 
depending on the magnitude of a. Note that very large values of a will affect the CFL 
stability condition of any explicit method regardless of the use of local time-stepping. 
Nevertheless, the resulting scheme will be explicit, since D is essentially a diagonal 
matrix. Second, we can no longer take advantage of any inherent symmetry in time of 
the solution. Third, to avoid any loss of accuracy, we must carefully initialize the LTS 
scheme, which again is based on the highly efficient (two-step) leap-frog method. 
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The exact solution z(t) of (3.2 1 satisfies 



z(t + At)-2 z(t) + z(t - At) + y D (z(t + At) - z{t - At)) 
= At 2 J (l-|0|)(R(t + 0At)-Az(t + 0Afc)) d6 + 0{At 4 ). 



(3.22) 



To derive a second-order LTS method for ( 3.2 1, we now split the vectors z(t) and R(t) 
as in (3.7 1 and approximate the integrands in (3.22 1 as follows: 

R [coarse] ^ _)_ Q At) + R [fine] (t + 6 At) ~ R[ coarse l (t) + PR(t + 6 At) , 

A (zl coarse ](t + At) + z[ flne l(i + 6 At)) ~ Az[ coarse l(t) + APz(6At) . 
We thus have 

z(t + At) - 2 z(t) +z{t- At) + y D (z(t + At)-z(t- At)) 
~ At 2 {(I - P)R(t) - A(I - P)z(t)} 

+ At 2 y (1 " {PR(i + #Ai) - APz(0Ai)} c/6» . 
Next for fixed t, let z(r) solve the differential equation 

0(r) + d|>) = (I - P)R(t) - A(I - P)z(t) + PR(i + r) 

- APz(r) , 
z(0) = z(i), z'(0) = v, 



(3.23) 



(3.24) 



where z/ will be specified below. Since the exact solution z(t) of (3.24 1 satisfies 



z(At) - 2 z(0) + z(-Ai) + y D (z(At) - z(-Ai)) 
= Ai 2 {(I - P)R(i) - A(I - P)z(t)} 

+ At 2 y (1 - \9\) {PR(t + 6 At) - APz(6At)} d6 , 



(3.25) 



from the comparison of ( |3.23[ l and p.25| , we infer that 

z(t + At) + z(t - At) + yD (z(t + At) - z(t - At)) 

~ z(At) + z(-Ai) + y D (z(At) - z(-At)) . 



(3.26) 
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In our local time-stepping scheme, we shall use the right side of (|3.26|l to approxi- 



mate the left side. In doing so, we must carefully choose v in (3.24 1 to minimize that 



approximation error. By using Taylor expansions and the fact that z and z solve (3.2 1 
and ( 3.24[ ), respectively, we obtain 



z(t + At) + z(t - At) = 2z(t) + z"(t)At 2 + 0{At 4 ) 

= 2z(t) + (R(t) - Az(t) - D z'(t))At 2 + 0(At 4 ) , 
z(At) +z(-At) = 2z(0) + z"(0)At 2 + 0(At 4 ) 

= 2z(t) + (R(t) - Az(t) - D v)At 2 + 0{At 4 ) , 

together with 

z(t + At) - z(t - At) = 2z'{t) At + 0(At 3 ) , z(At) - z(-At) =2vAt + 0{At 3 ) 



Hence for arbitrary v, the right side of ( 3.26 1 is not sufficiently accurate to approximate 
the left side while preserving overall second-order accuracy. However, if we choose 

v = z'(t) 



in ( 3.24 1, the 0(At ) terms in (3.26 1 cancel each other and overall second-order accu- 



racy of the scheme can be achieved. Since the term on the right side of (3.26 1 is not 



symmetric in time, unlike in the previous section (see (3.12 1 and (3.14l), we need to 
compute the value of z(r) both at r = At and at r = —At. 

For the numerical solution of ( 3.24| , we shall use the leap-frog scheme with the 
local time-step At = At /p. Since the leap-frog scheme is a two-step method, we need 
a second-order approximation of z'(0) = z'(t) during every initial local time-step. 
Since the value of z n+ i is still unknown at time t = t n , we now derive a second-order 
approximation z' n ~ z'(i) that uses only z„ and z„_i. First, we approximate 



J n-l/2 "n+1/2 



(3.27) 



where both z' 



n-l/2 



z'(t - At/2) and z' 



n+l/2 



z'(t + At/2) are second-order 



approximations. By using second-order central differences for z' 



s n-l/2 



Al 



and the differential equation (3.2 1 for z' r 



^n+l/2' 



ra-1/2' 



Zn-1 



+ 0(At 2 ) , 



(3.28) 



J n+l/2 ^n-1/2 

At 



+ Dz' r 



Az n + 0{At 2 ) , 
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we obtain 



I - ^Dj Z " J±± + AiR„ - AtAz„| + 0(At 2 ) . 

(3.29) 



Then, we insert (3.28 1, (3.29l into (3.27 1, which yields a second-order approximation 
of z'(t). 

In summary, the second-order LTS algorithm for the solution of ( |3.2[ ) computes 
z n+ i ~ z(t + At), for given z n and z„_i, as follows: 



LTS-LFCN2(p) Algorithm 

(i) Set w := (I — P)R„ — A(I — P)z n , zq := z n and 



1 



Zn-l 



A/ 



At 
I +y D 



I - -D ) Z " Z "" 1 + AtR„ - AtAz, 
2 At 



(ii) Compute 



'i/p '■= z o H z 



At , 1 /At\ 



2 VW 



(w + PR„o - APz - Dz' n ) awe? 



z_i/„ := z z 

if p 



At , 1 ( At\ 



(w + PR„, - APz - Dz'„ 

2 VP/ 



(hi) For m = 1, . . . ,p — 1, compute 



z (m+l)/p 



At 



2z 



m/p 



At\ 



2p l Z(m " 11 '" 



+ ( — J (w + PR„ m - APz m/p ) 



and 

z -(m+l)/p 



2p 



2z_ 



m/p 



At 

2p 



-(m-l)/p 



~J (w + PR„ _ m - APz_ m/p ) I . 



(iv) Compute 



At 
I + y D 



I-yD j i-z„ i +z ,) . 
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If a is piecewise constant in each element, M and M CT can be diagonalized si- 
multaneously and hence the matrix D is diagonal. If a varies in elements, D is a 
block-diagonal matrix and both (I ± (Ai/2p)D) and (I ± (Ai/2)D) can be explicitly 
inverted at low cost. In that sense, the LTS-LFCN2(p) scheme is truly explicit. Again, 
if the fraction of nonzero entries in P is small, the overall cost is dominated by the 
computation of w in step 1 . 

Proposition 3.4. For R(i) e C 2 ([0, T}), the local time-stepping method LTS-LFCN2 
is second-order accurate. 

Proof. See 123). □ 
Remark 3.5. For a = (D = 0), the LTS-LFCN2(p) algorithm coincides with the 



LTS-LF2(p) algorithm and thus also conserves the discrete energy (3.17l. For a ^ 
and p = 1, i.e. no local mesh refinement, one can easily show that the energy is 
no longer conserved but decays with time (independently of a) under the same CFL 
condition as in the case with a = 0. 



4 Adams-Bashforth based LTS methods for damped waves 



Starting from the standard leap-frog method, we proposed in Section[3]energy conserv- 
ing fully explicit LTS integrators of arbitrarily high accuracy for undamped waves. By 
blending the leap-frog and the Crank-Nicolson methods, a second-order LTS scheme 
was also derived there for damped waves, yet this approach cannot be readily extended 
beyond order two. To achieve arbitrarily high accuracy in the presence of damping, 
while remaining fully explicit, we shall derive here explicit LTS methods for damped 
wave equations based on Adams-Bashforth (AB) multi-step schemes. 



The H 1 -conforming and the IP-DG finite element discretizations of (2.1 1 presented 
in Section [2] lead to the second-order system of differential equations ( 2.4 », whereas 
the nodal DG discretization leads to the first-order system of differential equations 
( |2.1 1\ . In both (2.4 1 and ( |2.11| the mass matrix M is symmetric, positive definite and 
essentially diagonal; thus, M _1 or M - ? can be computed explicitly at a negligible 
cost. For simplicity, we restrict ourselves here to the homogeneous case, i.e. F(t) = 0. 

If we multiply (2.4 1 by M~5 ; we obtain (3.2 1. Thus, we can rewrite (3.2 1 as a 
first-order problem of the form 



dy 

dt 



(*) = By(t) , 



(4.1) 



with 



T 



B 
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Similarly, we can also rewrite (2.11 1 in the form (4.1 1 with y(t) = Q(t) and B = 
M _1 (—Mq- — C). Hence all three distinct finite element discretizations from Section 
[2] lead to a semi-discrete system as in (4.1 1. Starting from explicit multi-step AB 
methods, we shall now derive explicit LTS schemes of arbitrarily high accuracy for a 
general problem of the form ( |4.1| . 

First, we briefly recall the construction of the classical fc-step (fcth-order) Adams- 
Bashforth method for the numerical solution of (4.1 1 [29|. Let ti = iAt and y n , 
y n _i,..., y„_fc + i the numerical approximations to the exact solution y(t n ), 
y(£ n _fe + l). The solution of (4.1 1 satisfies 

r*n+£Ai 



y(£« + £Ai) = y(t n ) + / By(t)dt, < ^ < 1 



(4.2) 



We now replace the unknown solution y(t) under the integral in (4.2 1 by the interpo- 
lation polynomial p(t) through the points (ti, y{), i = n — k + 1, . . . , n. It is explicitly 
given in terms of backward differences 



vV 



V J+1 y„ = V J y„ - V^y 



n-l 



by 



fc-i 



P (t)=p(t n +sM) = Y,(-iy 

3=0 



Integration of (4.2 1 with y(t) replaced by p(t) then yields the approximation y n +£ of 
y{t n + £At),0<Z< 1, 



fc-i 



3=0 



(4.3) 



where the polynomials "fj($) are defined as 



7i(0 = 



ds . 



3 



They are given in Table [T] for j < 3. After expressing the backward differences in 
terms of y n -j an d setting £ = 1 in (4.3 1, we recover the common form of the fc-step 
Adams-Bashforth scheme [291 



fe-i 



y n +i =y n + AiB ajy 
i=o 



n-j , 



(4.4) 



where the coefficients otj, j — 0, . . . , k — 1 for the second, third- and fourth-order 
(fc = 2, 3,4) Adams-Bashforth schemes are given in Table [2] For higher values of fc 
we refer to J29). 
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3 


12 3 


7i(0 


t 1 c2 lrf i U2 1 <4 , 1 13 i 1 tZ 
? 2? 6? + 4? 24? + 6? + 6? 



Table 1. Coefficients 7j(£) for the explicit Adams-Bashforth methods. 





a 


a i 


Q 2 


°3 


k 


= 2 


3 
2 


1 

2 








k 


= 3 


23 
12 


16 
12 


5 

12 





k 


= 4 


55 
24 


59 
24 


37 
24 


9 

24 



Table 2. Coefficients for the A:-th order Adams-Bashforth methods. 



Starting from the classical AB methods, we shall now derive LTS schemes of arbi- 
trarily high accuracy for ( |4.1[ >, which allow arbitrarily small time-steps precisely where 
small elements in the spatial mesh are located. To do so, we first split the unknown 
vector y(t) in two parts 

y(t) = (I - P)y(t) + Py(t) = y^(t) + y^(t) , 

where the matrix P is diagonal. Its diagonal entries, equal to zero or one, identify 
the unknowns associated with the locally refined region, where smaller time-steps are 
needed. 

The exact solution of ( |4.1[ l again satisfies 

y(t n + £At) = y(t n ) + / B (y^(t) + y^(t)) dt , < £ < 1 . 

(4.5) 

Since we wish to use the standar d fc-s tep Adams-Bashforth metho d in the coarse re- 
gion, we approximate the term in (4.5 1 that involve y[ coaise l (t) as in (4.2 1, which yields 



y(t„ + fAt) « y n + At B(I - P) 7i(0V J yn + / Bp y W dt • ( 4 -°) 

j=0 *" 



To circumvent the severe stability constraint due to the smallest elements associated 
with y [finel (i), we shall now treat y [flne] (i) differently from y[ coarse l(i). Hence, we 
instead approximate the integrand in (|4.6|l as 



rt n +£At r£At 

/ BPy(t)dt« / BPy(r)dr, 

Jt n Jo 
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3 


12 3 


7,(0 





Table 3. The polynomial coefficients 7j(£) 



where y(r) solves the differential equation 



k-l 



<,y 'r) = B(I - P) ^7, (^) V'y„ + BPy(r) , 



dr y ' y ' ^ u \At 

3=0 



(4.7) 



y(0) = y n , 

with coefficients 



7i(f) = ^7,(0 = (W f f da) = (~iy( f J . (4.8) 



The polynomials 7j(£) are given in Table|3]for j < 3. Replacing y(t) by y(t) in (4.6 1, 
we obtain 



fc-i 



y(t„ + £At)«y„+AtB(I-P) V>i(£)V J y„ + / BPy(r)dr. (4.9) 
By considering ( |4.7| l in integrated form, we find that 



y(£Ai)=y(0) + B(I-P)£ / ^ (-J dr V J y„ + / BPyMdr 

= y J1 + AtB(I-P)V7 J (0V J y„+ / BPy(r)rfr. 

^0 



(4.10) 

From the comparison of (4.9 1 and ( |4.10| we infer that 

y(t n + £At) « y(£At) . 



Thus to advance y(t„) from t n to t„ + At, we shall evaluate y (At) by solving (4.7 1 on 
[0, At] numerically. We solve (4.7 1 until r = At again with a /c-step Adams-Bashforth 
scheme, using a smaller time-step At = At/p, where p denotes the ratio of local 
refinement. For m = 0, . . . ,p — lwe then have 

k ~ l k ~ l „ fm-l\ 
+ At B(I - P) J] a € ]T 7, ( — ) V'y 

£=0 j=0 

fc-i 

+ ArBP^a £ y (m _; )/! ,, 
£=0 



P 

(4.11) 
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where ag, £ = 0, . . . , k — 1 denote the coefficients of the classical fc-step AB scheme 
(see Table[2]l. Finally, after expressing the backward differences in terms of y n -e, we 
find 

fe-l k-l 

y(m+i)/p = y m / P + At B(I - P) Pm,i y n -e + At BP ^ a g y (m _ z) / p , (4.12) 
where the constant coefficients (3 m ^, m — 0, . . . ,p — 1, £ — 0, . . . , k — 1, satisfy 

/W = E^D- 1 )'(j) ; w( !! y i ). (4-13) 



with 7j defined in (4.8 1. 



j ... j 



In summary, the LTS-AB/c(p) algorithm computes y n +i — y(t n + At), given y n 

Yn— ls*'*s yn— /e+l» 

B(I-P)y n _ 1 ,...,B(I-P) 

y?i— fc+i 

and Py„_i/ P , Py„_ 2 / P 

Py„_ (A; _ 1) / p as follows: 

LTS-AB/c(p) Algorithm 

(i) S<?r y := y n , Y-i/ P ■= Py n -e/ P > £=!,-■■, k-l. 

(ii) Setw n _ e :=B(I-P)y n _^ = l,...,fc-1. 
(hi) Compute w n := B(I — P)y n . 
(iv) For m — 0, . . . ,p — 1, compute 



., fe-i A . fc-i 

At x-^ „ At, 



y(m+i)/ P : = y m / P h — 2^ ^ m ' £ w "-^ — af y ( 



p tf P 1=0 



m—l)/p ■ 



(v) Sety n+l := y h 

Steps 1-4 correspond to the numerical solution of ( |4.7| until t = At with the fe- 
step AB scheme, using the local time-step At = At/p. For P = or p = 1, that is 
without any local time-stepping, we thus recover the standard fc-step Adams-Bashforfh 
scheme. If the fraction of nonzero entries in P is small, the overall cost is dominated 
by the computation of w„ in Step 3, which requires one multiplications by B(I — P) 
per time-step At. All further matrix-vector multiplications by BP only affect those 
unknowns that lie inside the refined region, or immediately next to it; hence, their 
computational cost remains negligible as long as the locally refined region contains a 
small part of £2. 

We have shown above how to derive LTS-AB/c(p) schemes of arbitrarily high ac- 
curacy. Since the third- and fourth-order LTS-AB/c(p) schemes are probably the most 
relevant for applications, we now describe the US-ABk(p) schemes for k = 3, 4 and 
p = 2. Other examples of LTS Adams-Bashforfh schemes are listed in ll26ll . 
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For k = 3 and p = 2, the LTS-AB3(2) method reads: 



yi/2 



yn+i = yi = yi/2 



A/ 

y 

A/ 

y 

A/ 

y 

A/ 



B(I-P) 

"23 
BP — y 

B(I-P) 

~23„ 



17 7 2 

j^y« — j^yn— i t Y^y«— 2 



16~ 5 

Y^y-1/2 + ^y™- 1 



29 25 8 

Yn — T^ra-l + 7^yra-2 



12' 



12' 



12" 



BP 



16 



Y^-yi/2 ~~ i2^ n ~*~ 12 / 2 



For the case with fc = 4 and p = 2, we find the LTS-AB4(2) scheme: 



yi/2 



A/. 



B(I-P) 



+ y BP 

At, 



55 
24" 



297 
19T 

59. 



187 



107 



25 



192 y " _1 + 192 y "" 2 192 y ™" 3 



37 



24 -y-i/ 2 + ^yn-i 



24 



y-3/2 



y«+i = yi = yi/2 + y B(l - P) 



583 



757 



485 



119 



192 y " 192 y " 1 192 y ™ 2 192 y " 3 



A* „ 
+ y BP 



55. 



59 



37. 



24 yi / 2 ~ 24 y ™ 24 y_1/2 ~ 24 y ™ _1 



Proposition 4.1. The local time-stepping method LTS-ABk(p) is consistent of order k. 
Proof. See G6l. □ 



5 Numerical results 

Here we present numerical experiments that validate the expected order of conver- 
gence of the above LTS methods and demonstrate their usefulness in the presence of 
complex geometry. First, we consider a simple one-dimensional test problem illus- 
trate the stability properties of the different LTS schemes presented above and to show 
that they yield the expected overall rate of convergence when combined with a spatial 
finite element discretization of comparable accuracy, independently of the number of 
local time-steps p used in the fine region. Then, we illustrate the versatility of our LTS 
schemes by simulating the propagation of a circular wave in a square cavity with a 
small sigma- shaped hole. 



5.1 Stability 

We consider the one-dimensional homogeneous damped wave equation <\2.l\ with con- 
stant wave speed c = 1 and damping coefficient a = on the interval D, = [0,6]. 
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Next, we divide D, into three equal parts. The left and right intervals, [0 , 2] and [4 , 6], 
respectively, are discretized with an equidistant mesh of size /i coaise , whereas the in- 
terval D,f = [2 , 4] is discretized with an equidistant mesh of size h Rne = /i coarse /p. 
Hence, the two outer intervals correspond to the coarse region whereas the inner inter- 
val [2 , 4] to the refined region. In lPT6il . we have studied numerically the stability of the 
LTS-LF2(p) and the LTS-LFME4(p) methods. To determine the range of values At for 
which the LTS-LF2(p) scheme is stable, the eigenvalues of (Ai 2 /4)A P (A p is defined 
by ( 3.16 )) for varying At/Atip are computed, where Atpp denotes the largest time- 
step allowed by the standard leap-frog method. The LTS-LF2(p) scheme is stable for 
any particular At if all corresponding eigenvalues lie between zero and one; otherwise, 
it is unstable. We have observed that the largest time step allowed by the LTS-LF2(p) 
scheme is only about 60% of Atpp. A slight extension (overlap) of the region where 
local time steps are used into that part of the mesh immediately adjacent to the refined 
region typically improves the stability of the LTS-LF2(p) scheme. Moreover, the nu- 
merical results suggested that an overlap by one element when combined with a "P 1 
continuous FE discretization (with mass lumping), or by two elements when combined 
with a IP-DG discretization, permits the use of the maximal (optimal) time step Atpp. 
The numerical results also suggested that an overlap by one element for the IP-DG 
discretization is needed for the optimal CFL stability condition of the LTS-LFME4(p) 
independently of p. Remarkably, no overlap is needed for the LTS-LFME4(p) scheme 
to remain stable with the optimal time- step when combined with the continuous V 3 
elements. 

In ll26l . we have considered the above one-dimensional problem with a = 0.1. 
We have written the LTS-ABfc(p) scheme as a one-step method and than studied nu- 
merically it stability when combined with a spatial finite element discretization of 
comparable accuracy. For a spatial discretization with standard continuous, IP-DG or 
nodal DG finite elements, we have obtained that the maximal time-step At p allowed 
by the LTS-AB 2(p) scheme is about 80 % of the optimal time-step AtABi (the largest 
time-step allowed by the standard two-step AB method) independently of h, p and cr; 
moreover, the CFL stability condition of the LTS-AB 3{p) and LTS-AB 4(p) schemes 
is optimal for all h, p and a. 



5.2 Convergence 



We consider the one-dimensional homogeneous model problems (2.1 1 and (2.9 1 with 
constant wave speed c = 1 and damping coefficient a = 0.1 on the interval D, = 
(0 , 6). The initial conditions are chosen to yield the exact solution 



u(x, t) 



2e~~ 



V4n 2 - a 2 
du 



. sin 



(irx) sin ( - \/ 4-ir 2 — a 2 ] , 



v(x,t) = —(x,t) , w(x, t) = - Vu{x, t) . 
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(a) continuous FE (h = 0.2, 0.1, 0.05, 0.025) (b) IP-DG (h = 0.2, 0.1, 0.05, 0.025) 




(c) nodal DG (h = 0.02, 0.01, 0.005, 0.0025) 



Figure 1. LTS-AB4(p) error vs. h = ft coarse for V 3 finite elements with p = 2, 5, 7. 



Again, we divide D, into three equal parts. The left and right intervals, [0,2] and 
[4, 6], respectively, are discretized with an equidistant mesh of size ft coarse , whereas the 
interval [2,4] is discretized with an equidistant mesh of size /i flne = /i coarse /p. Hence, 
the two outer intervals correspond to the coarse region and the inner interval [2 , 4] to 
the refined region. 

First, we consider a V 3 continuous FE discretization with mass lumping and a se- 
quence of increasingly finer meshes. For every time-step At, we shall take p > 2 local 
steps of size At = At/p in the refined region, with the fourth-order time-stepping 
scheme LTS-AB4(p). The first three time-steps of each LTS-AB4(p) scheme are ini- 
tialized by using the exact solution. According to our results on stability, we set 
At = AtAB4, the corresponding largest possible time-step allowed by the AB ap- 
proach of order four on an equidistant mesh with h = /i coarse . As we systematically 
reduce the global mesh size /i coarse , while simultaneously reducing At, we monitor the 
L 2 space-time error in the numerical solution \\u(-,T) — u h (-,T)\\ L 2^ at the final 
time T = 10. In frame (a) of Fig. [T] the numerical error is shown vs. the mesh size 
h = /! coarse ; regardless of the number of local time-steps p = 2, 5 or 7, the numerical 
method converges with order four. 

We now repeat the same experiment with the IP-DG (a = 20 in (|2.7[l) and the nodal 
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Figure 2. Two-dimensional example: the computational domain Q,. 



DG discretizations with "P 3 -elements. As shown in frames (b) and (c) of Fig. [I] the 
LTS-AB4(p) method again yields overall fourth-order convergence independently of 
P- 

Remark 5.1. We have obtained similar convergence results for other values of p and a. 
In summary, we observe the optimal rates convergence of order k for the LTS-AB/c(p) 
schemes as well as for the LTS-LF2(p), LTS-LSME4(p) and LTS-LFCN2(p) schemes, 
regardless of the spatial FE discretization and independently of the number of local 
time-steps p and the damping coefficient a. For more details, we refer to iTR)! 1251127 



5.3 Two-dimensional example 



To illustrate the usefulness of the LTS method presented above, we consider ( 2. 1 1 in a 
square cavity D. = (0, l) 2 , with a small sigma-shaped hole - see Figure [5j We set the 
constant wave speed c = 1 and the damping coefficient 

/ 10, x 2 <0.5 
0.1 , otherwise . 



We impose homogeneous Neumann conditions on the boundary of D, and choose as 
initial conditions 



u (x) 



exp 



xo|| 2 A 2 ) 



u (x) = 0. 



where x = (0.45,0.55) and r = 0.012. 

For the spatial discretization we opt for the V 2 continuous finite elements with mass 
lumping. First, D, is discretized with triangles of minimal size /j, coarse = 0.03. How- 
ever, such triangles do not resolve the small geometric features of the sigma-shaped 
hole, which require 

h hne _ ^coarse / 7j as shown j n Figure [5] Then, we successively 
refine the entire mesh three times, each time splitting every triangle into four. Since 
the initial mesh in D, is unstructured, the boundary between the fine and coarse mesh 
is not well-defined. Given /i coarse ; here the fine mesh corresponds to all triangles with 
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Figure 3. The triangular initial mesh at various magnification rates: the darker triangles 
belong to the "fine" mesh. 




Figure 4. Gaussian pulse penetrating a cavity with a small sigma-shaped hole. The 
solution is shown at times t = 0.1, 0.2, 0.3, 0.4, 0.44 and 0.5. 

h < 0.75/i coarse in size, that is the darker triangles in Figure [5] The corresponding 
degrees of freedom in the finite element solution are then selected merely by setting to 
one the corresponding diagonal entries of the matrix P. 

For the time discretization, we choose the third-order LTS-AB3(7) time-stepping 
scheme with p — 7, which for every time-step At takes seven local time-steps inside 
the refined region. Thus, the numerical method is third-order accurate in both space 
and time under the CFL condition At = 0.07 /i coarse ? determined experimentally. If 
instead the same (global) time-step At was used everywhere inside £2, it would have to 
be about seven times smaller than necessary in most of D,. As a starting procedure, we 
employ a standard fourth-order Runge-Kutta scheme. 

In Fig. |4j snapshots of the numerical solution are shown at different times. The 
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circular wave, initiated by the Gaussian pulse, propagates outward until it impinges 
first on the sigma-shaped hole and later on the upper and left boundaries of CI. The 
reflected waves move back into D, while multiple reflections occur both at the obstacle 
and along the interface at X2 = 0.5. As the waves cross that interface and penetrate 
the lower part of D., they are strongly damped. 

6 Concluding remarks 

Starting from the classical leap-frog (LF) or Adams-Bashforfh (AB) methods, we have 
presented explicit local time-stepping (LTS) schemes for wave equations, either with 
or without damping. By allowing arbitrarily small time-steps precisely where the 
smallest elements in the mesh are located, these LTS schemes circumvent the crip- 
pling effect of locally refined meshes on explicit time integration. 

When combined with a spatial finite element discretization with an essentially diag- 
onal mass matrix, the resulting LTS schemes remain fully explicit. Here three such fi- 
nite element discretizations were considered: standard H l -conforming finite elements 
with mass-lumping, an IP-DG formulation, and nodal DG finite elements. In all cases, 
our numerical results demonstrate that the resulting fully discrete numerical schemes 
yield the expected space-time optimal convergence rates. Moreover, the LTS-AB(/c) 
schemes of order k > 3 have optimal CFL stability properties regardless of the mesh 
size h, the global to local step-size ratio p, or the dissipation a. Otherwise, the CFL 
condition of the LTS scheme may be sub-optimal; then, by including a small overlap 
of the fine and the coarse region, the CFL condition of the resulting LTS scheme can 
be significantly enhanced. 

Since the LTS methods presented here are truly explicit, their parallel implemen- 
tation is straightforward. Let At denote the time-step imposed by the CFL condition 
in the coarser part of the mesh. Then, during every (global) time-step At, each local 
time-step of size At/p inside the fine region of the mesh simply corresponds to sparse 
matrix-vector multiplications that only involve degrees of freedom associated with the 
fine region of the mesh. Those "fine" degrees of freedom can be selected individually 
and without any restriction by setting the corresponding entries in the diagonal pro- 
jection matrix P to one; in particular, no adjacency or coherence in the numbering of 
the degrees of freedom is assumed. Hence the implementation is straightforward and 
requires no special data structures. 

In the presence of multi-level mesh refinement, each local time-step in the fine re- 
gion can itself include further local time-steps inside a smaller subregion with an even 
higher degree of local mesh refinement. The explicit local time-stepping schemes de- 
veloped here for the scalar damped wave equation immediately apply to other damped 
wave equations, such as in electromagnetics or elasticity; in fact, they can be used for 
general linear first-order hyperbolic systems. 

Acknowledgments. We thank Manuela Utzinger and Max Rietmann for their help 
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